home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / polar_surface.pro < prev    next >
Text File  |  1997-07-08  |  3KB  |  98 lines
  1. ; $Id: polar_surface.pro,v 1.5 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1992-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5.  
  6. function polar_surface, z, r, theta, GRID = grid, SPACING = sp, $
  7.     BOUNDS = bounds, QUINTIC = quintic, MISSING = missing
  8. ;+
  9. ; NAME:
  10. ;    POLAR_SURFACE
  11. ;
  12. ; PURPOSE:
  13. ;    This function interpolates a surface from polar coordinates
  14. ;    (R, Theta, Z) to rectangular coordinates (X, Y, Z).
  15. ;
  16. ; CATEGORY:
  17. ;    Gridding.
  18. ;
  19. ; CALLING SEQUENCE:
  20. ;    Result = POLAR_SURFACE(Z, R, Theta)
  21. ;
  22. ; INPUTS:
  23. ;    Z:     An array containing the surface value at each point.
  24. ;         If the data are regularly gridded (GRID=1) in R and 
  25. ;         Theta, Z is a two dimensional array, where Z[i,j] has a
  26. ;         radius of R[i] and an azimuth of Theta[j].  If GRID is
  27. ;         not set, R[i] and Theta[i] contain the radius and azimuth
  28. ;         of each Z[i].
  29. ;    R:     The radius. If GRID is set, Z[i,j] has a radius of R[i].
  30. ;         If GRID is not set, R must have the same number of elements
  31. ;         as Z, and contains the radius of each point.
  32. ;    Theta:   The azimuth, in radians. If GRID is set, Z[i,j] has an
  33. ;         azimuth of Theta[j]. If GRID is not set, Theta must
  34. ;         have the same number of elements as Z, and contains
  35. ;         the azimuth of each point.
  36. ;
  37. ; KEYWORD PARAMETERS:
  38. ;    GRID:    Set GRID to indicate that Z is regularly gridded in R
  39. ;         and Theta.
  40. ;    SPACING: A two element vector containing the desired grid spacing
  41. ;         of the resulting array in X and Y.  If omitted, the grid
  42. ;         will be approximately 51 by 51.
  43. ;    BOUNDS:  A four element vector, [X0, Y0, X1, Y1], containing the
  44. ;         limits of the XY grid of the resulting array.  If omitted,
  45. ;         the extent of input data sets the limits of the grid.
  46. ;    QUINTIC: If set, the function uses quintic interpolation, which is
  47. ;         slower but smoother than the default linear interpolation.
  48. ;    MISSING: A value to use for areas within the grid but not within
  49. ;         the convex hull of the data points. The default is 0.0.
  50. ;
  51. ; OUTPUTS:
  52. ;    This function returns a two-dimensional array of the same type as Z.
  53. ;
  54. ; PROCEDURE:
  55. ;    First, each data point is transformed to (X, Y, Z). Then
  56. ;    the TRIANGULATE and TRIGRID procedures are used to interpolate
  57. ;    the surface over the rectangular grid.
  58. ;
  59. ; EXAMPLE:
  60. ;    r = findgen(50) / 50.                  ;Radius
  61. ;    theta = findgen(50) * (2 * !pi / 50.)         ;Theta
  62. ;    z = r # sin(theta)        ;Make a function (tilted circle)
  63. ;    SURFACE, POLAR_SURFACE(z, r, theta, /GRID)      ;Show it
  64. ;
  65. ; MODIFICATION HISTORY:
  66. ;    DMS     Sept, 1992    Written
  67. ;-
  68.  
  69.  
  70. IF keyword_set(grid) THEN BEGIN        ;Regulary gridded?
  71.     s = size(z)
  72.     if s[0] ne 2 then message, "Z must be 2d if GRID is set"
  73.     if n_elements(r) ne s[1] then $
  74.         message, "R has wrong number of elements"
  75.     if n_elements(theta) ne s[2] then $
  76.         message, "Theta has wrong number of elements"
  77.     x = r # cos(theta)
  78.     y = r # sin(theta)
  79. ENDIF ELSE BEGIN
  80.     if (n_elements(r) ne n_elements(z)) or $
  81.         (n_elements(theta) ne n_elements(z)) then $
  82.         message, "R and Theta must have as many elements as Z"
  83.     x = r * cos(theta)
  84.     y = r * sin(theta)
  85. ENDELSE
  86.  
  87. TRIANGULATE, x, y, tr        ;Triangulate it
  88. if n_elements(bounds) lt 4 then $
  89.     bounds = [ min(x, max=x1), min(y, max=y1), x1, y1]
  90. if n_elements(sp) lt 2 then $
  91.     sp = [ bounds[2] - bounds[0], bounds[3] - bounds[1]] / 50.
  92. if n_elements(missing) le 0 then missing = 0
  93.  
  94. RETURN, TRIGRID(x, y, z, tr, sp, bounds, QUINTIC = KEYWORD_SET(quintic), $
  95.     MISSING = missing)
  96. END
  97.  
  98.